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SUMMARY 


Mathematical development and some computed results are presented for 
Mindlin plate and shell elements, suitable for analysis of laminated composite 
and sandwich structures. These elements use the conventional 3 (plate) or 5 
(shell) nodal degrees of freedom, have no communicable mechanisms, have 
no spurious shear energy (no shear locking), have no spurious membrane 
energy (no membrane locking) and do not require arbitrary reduction of out- 
of-plane shear moduli or under-integration. Artificial out-of -plane rotational 
stiffnesses are added at the element level to avoid convergence problems or 
singularity due to flat spots in shells. 

In regular rectangular meshes, the Martin-Breiner 6-node triangular curved 
shell (MB6) is about equivalent to the conventional 8-node quadrilateral with 
2x2 integration (which is quite good). In distorted meshes, this 6-node 
triangular element is distinctly better. The accuracy of the MB6 is most 
evident in the NAFEMS LE2 curved shell patch test, where error at the 
specified point is only 0. 12 percent, and maximum error anywhere in the 
patch is 2.5 percent. In contrast, results for five 6 and 8 node elements in 
commercial programs showed errors of 8 percent to 85 percent at the 
specified point, and maximum errors as large as -99 to +100 percent at some 
points in the patch. The four-node quadrilateral, MB4, has very good 
accuracy for a four-node element, and may be preferred in vibration analysis 
because of narrower bandwidth. 

The mathematical developments used in these elements, included here in 
seven appendices, have been applied to elements with 3, 4, 6, and 10 nodes 
and can be applied to other nodal configurations. 


INTRODUCTION 

Since the inception of finite element analysis, efforts have been made to 
develop accurate shell elements. Many formulations have been tried, and no 
attempt at review will be made here. Some element formulations are plagued 
by spurious shear strains (“shear locking”); some by spurious membrane 
strains (“membrane locking”). Some elements or element options may 
represent bending deformation well but not membrane deformation. Other 
elements or element options may represent membrane deformation well but 
not bending. Few have been equally accurate for all deformation modes. The 



unwary analyst may assume that because a given element and mesh solves 
one load case correctly, it will be equally accurate for a different load case. 
The eight- node isoparametric shell with 2x2 integration may be quite 
accurate when rectangular, but it is under-integrated and much less accurate 
when its shape is distorted. 

The six-node curved triangular shell, MB6, presented here is believed to be a 
significant advancement because it has excellent and uniform accuracy in all 
deformation modes, and needs no “options” for different load cases, etc. The 
triangular shape allows easy mesh generation around openings and 
discontinuities where a rectangular element’s shape must be severely 
distorted. Accuracy of stresses is most improved where it is most important; 
in regions of high stress gradients. 

The four-node warped quadrilateral shell, MB4, performs very well, and may 
be preferred for vibration analysis because of its narrower bandwidth. In- 
plane deformation can be improved by optional inclusion of incompatible 
modes [4]. Bending deformation can be improved by optional activation of a 
simulated antisymmetric bending mode. 


ELEMENT PERFORMANCE 

The National Agency for Finite Element Methods and Standards in the U.K. 
supported development of testing procedures and test cases for evaluating 
finite elements and programs, beginning with a set of benchmark problems 
by Barlow and Davies in 1986 [5]. Four of these NAFEMS test cases are 
shown in Figures 1,2,6, and 7. In addition, some other popular test cases are 
shown; the Scordelis-Lo roof [1,2] in Figure 3; the Morley skew plate [6] in 
Figure 4; and Plunkett’s vibrating wedge [7] in Figure 5. Some results for 
these test cases are shown in Tables 1 to 7. 

Results of NAFEMS tests for some commercial programs were published in 
Benchmark Magazine in 1989 [8]. Some results for the LE2 and LE3 shell 
test cases are shown in Table 8. Errors were quite significant, and at least one 
entry appears over-optimistic. 

More recent calculations (February 1999) suggest that some commercial 
programs still contained shell elements with poor accuracy. Table 9 shows 
results for 5 important test cases for six-node elements in two versions of 
NASTRAN, and for the Ahmad-Irons-Zienkiewicz (AIZ) [1] elements which 
are contained in some popular programs, as well as the MB6. The NAFEMS 
LE2 test, described in Figure 1, only asks for the stress at the outer surface at 
point E, and errors in these average stresses are shown in Table 9. Some 
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stresses at other points in the patch are shown in parentheses. Note that this is 
a patch test, and stress should be constant everywhere. However, it is possible 
to get stresses from 99 percent low to about 100 percent high from one of 
these commercial programs at some points in the patch. Selection of the 
“best” element option may improve results, but criteria for selecting these 
options and reasons for which is the default are not necessarily clear. In 
contrast, the MB6 error in the LE2 test is only 0.12 percent at the specified 
point, and its maximum error anywhere in the patch for either load case is 2.5 
percent. In addition, the default element option in “NASTRAN B” did not 
converge for either 4 or 6 node elements in the Scordelis-Lo roof [1,2], 

Figure 3. 

The MB 6 and MB4 elements perform very well on the critical test cases in 
Figures 1 to 7. 

Some results are presented in tables 1 to 7 from seven of the test cases 
used in element validation. Elements are identified as follows: 

MB 10 is the Martin-Breiner 10-node triangular shell 

MB6 is the Martin-Breiner 6-node triangular shell 

MB4 is the Martin-Breiner 4-node quadrilateral shell 

MB4SI is the Martin-Breiner 4-node quadrilateral with simulated 

antisymmetric bending mode and incompatible in-plane modes active. 

AIZ6 is the Ahmad-Irons-Zienkiewicz 6-node triangular shell. 

AIZ8 is the Ahmad-Irons-Zienkiewicz 8-node quadrilateral shell. 

AIZ10 is the Ahmad-Irons-Zienkiewicz 10-node triangular shell 
LU71S is a 3-node triangular shell [3, 10, 16] with the addition of a 
simulated antisymmetric bending mode (SABM). 

In NAFEMS LE2 [5], Figure 1 and Table 1, the stress should be uniform 
throughout. MB 6 is the only element tabulated which has a near-zero error in 
average stress and near-zero standard deviation in both cases. The MB6 
stress is also nearly uniform throughout, whereas for some other elements it 
is not. The maximum error for the MB6 at any node of any element, either 
top or bottom surface and either load case is 2.5 percent. 

In NAFEMS LE3 [5], Figure 2 and Table 2, the MB6 is much more accurate 
than the other quadratic elements in the coarse mesh. The LU71S element is 
surprisingly accurate in the coarse mesh. 

In the Scordelis-Lo Roof [1,2], Figure 3 and Table 3, the MB6 has, by far, the 
fastest convergence. In comparing calculated results, divergence was 
observed in some other elements/programs at about 16 nodes per side. 
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In the Morley skew plate, [6], Figure 4 and Table 4, the MB6 has very good 
accuracy, and other elements shown require many more degrees of freedom 
to get under 1 percent error. 

Plunkett’s vibrating wedge [7], Figure 5 and Table 5, is a severe test of 
performance with variable thickness. The MB6 computed the most (10) 
mode shapes that corresponded to Plunkett’s sketches, whereas the AIZ8 
computed only 8. Also, the MB6 errors in frequencies were smaller than the 
AIZ8 in 5 of those first 8 modes, and smaller than the AIZ6 in 9 of the 10 
modes. 

In NAFEMS T1 [5], Figure 6 and Table 6, both 6-node shells MB6 and AIZ6 
give excellent results, whereas the 8-node shell, AIZ8, is 7.8 percent low. 
Performance of the 4-node shell, MB4, with incompatible in-plane modes is 
very good. 

In the NAFEMS laminated strip test case, [9] Figure 7 and Table 7, all of the 
elements compared give very accurate results for deflection and bending 
stress. The MB6 error in interlaminar stress is 2.4 percent, which is good, 
although the AIZ element errors are even smaller. 

It is expected that the 6-node triangle, MB6, will be a particularly useful 
element in stress analysis of laminated composite plate and shell structures. It 
has nearly the same excellent accuracy in all deformation modes, and needs 
no “options” for different conditions. It needs only 3 integration points, and 
thus can be computationally efficient. Two triangles could easily be joined to 
make a 9-node or 8-node quadrilateral with a total of 6 integration points. 
Triangular elements have obvious advantages in modeling around 
discontinuities such as openings, joints and reinforcements where stresses are 
highest and most important. 


MATHEMATICAL DEVELOPMENT 


The basic formulation of the elements described here follows Ahmad-Irons- 
Zienkiewicz [1] with some modifications. The basic formulation is reviewed 
in Appendix A. Note that matrices are generated in element coordinates, with 
nodes 1,2,3 defining the x-y plane. One significant change is that nodal 
rotations are interpolated as expressed in Equation A5. It follows that 
bending strain terms do not include terms with derivatives of thickness, as in 
Equation A 12, since they would produce strains due to rigid body motion. 
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This change was made to avoid singularity or near singularity in tapered 
elements with very thin edges, as in Plunkett’s vibrating wedge [7]. If only 
constant thickness elements or elements with modest taper are to be analyzed, 
this modification is not necessary or even helpful. All of the numerical results 
presented here were generated by computer code that builds the element 
matrix layer-by-layer with two integration points through the thickness of 
each layer. However, it is shown in Equations A17 to A19 that code could be 
written with explicit integration through the thickness, which should execute 
faster. Note that only the stretching and bending parts of the stiffness are 
computed by the process in Appendix A, and an important modification to 
the membrane strain is described in Appendix B. The out-of-plane shearing 
stiffness is computed by a different process, shown in Appendix C. 


Appendix B shows an important innovation, not previously published, which 
eliminates spurious membrane strain or “membrane locking”. In curved 
quadratic elements (e.g. the AIZ 6-node triangle and 8-node quadrilateral) 
bending causes spurious mid-surface strains except at the 2x2 Gauss points. 
These spurious strains cause serious errors in fully integrated AIZ elements 
when the “rise” (deviation from flatness) of the element is only about 1/5 of 
the thickness. A happy exception is the (under-integrated) 8-node rectangular 
quadrilateral with 2x2 integration, since the spurious strains are correctly zero 
at the 2x2 Gauss points. However, this is little help in a triangular element or 
in a quadrilateral whose shape is significantly distorted. 

The technique used here to eliminate spurious membrane strain is, in concept, 
surprisingly simple. It is observed that the average mid surface strain due to 
constant-moment bending in the element is correct. The correct average strain 
can be obtained by averaging strains at the integration points. Unfortunately, 
this eliminates the gradient in mid-surface strains, which is needed for 
accurate solution of some problems, and it introduces mechanisms. 

In the 6-node triangular element, the gradient can be restored by using strains 
from triangular sub-regions, as shown in Appendix B. It can be shown that 
this process for recovering the mid surface strain gradient is exact for all 
constant strain states and all linear strain states in a flat 6-node triangle with 
straight sides. 

This concept probably could be applied to the 8-node quadrilateral with full 
(3x3) integration, which should make it capable of more distortion in shape. 
However, a 9-node or 8-node quadrilateral with only 6 total integration points 
can easily be generated by joining two 6-node triangles, which should be 
about equally accurate and require less compute time. 


Appendix C shows the method for calculating out-of-plane shear strain. This 
concept originated with Utku [3], who applied it to a 3-node flat shell. The 
basic concept is that the function used to interpolate out-of-plane 
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displacement must be one order higher than the function used to interpolate 
rotations. This eliminates the problem of spurious shear strains (“shear 
locking”). The concept was applied to laminated composite shells in 
Reference [10], extended to a 6-node triangle by Yu [1 1], and generalized to 
any nodal configuration by Martin and Breiner [12]. 

The procedure in Appendix C is strictly valid only for flat elements. 
However, the MB6 matrices are generated in element axis with the x-y plane 
defined by nodes 1-2-3. If the included angle of the element (angle between 
outward normals at opposite sides) is 20 degrees, the angle between the 
outward normal and the z-axis does not exceed 10 degrees. The cosine of 10 
degrees is 0.985, so the error is small. This could be considered a limitation 
or the MB6; that its included angle should not exceed about 20 degrees. 
However, in the Scordelis-Lo roof [1,2], Figure 3 and Table 3, the MB 6 
solution is quite accurate with only two elements, which span 40 degrees. 
More study of this limitation is needed, and the method should be developed 
for highly curved elements. 

Elements that use this procedure have one mechanism, which however is 
suppressed by joining two elements. The mechanism can be physically 
described as a relative rotation of top and bottom surfaces about the centroid, 
with no strain energy. Although deflections are always quite accurate, some 
imperfections in element strains in the MB 6 have been observed in a linear 
bending problem. Accuracy of strains is improved by using the weighted 
least squares fit of equation C20 andC21. This may be related to the 
mechanism, and more study of this may be appropriate. 


Appendix D shows the method of calculating strains at nodes used in the new 
elements. 


Appendix E shows the method, not previously published, used in all elements 
to generate artificial rotational stiffness about the z-axis. This is necessary 
because these elements inherently have only two rotational stiffness degrees 
of freedom at each node. The artificial stiffness must be added so that 
coordinate transformation of the element matrix and assembly with 6 degrees 
of freedom per node is possible. Selecting these artificial stiffnesses so that 
they avoid mechanisms, constraints and ill conditioning, but do not 
significantly stiffen the structure or add strains due to rigid body motion has 
been a persistent problem. 

The method discussed here meets all of the criteria just mentioned. It is a 
combination of one reported by Zienkiewicz [13] and one due to Kanok 
Nukulchai [14]. In a flat plate, each of these methods requires fixing at least 
one out-of-plane rotational degree of freedom to avoid a singularity. Used in 
combination, there are no mechanisms, and no special attention is needed for 
calculations such as vibration of a flat free-free plate. 
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Appendix F shows the development of the simulated antisymmetric bending 
mode (SABM) which is an option in 3 and 4-node elements. Linear Mindlin 
elements can express symmetric bending exactly, but are not capable of 
antisymmetric bending. The SABM substitutes shear deformation for 
antisymmetric bending deflection with the objective of preserving the correct 
total strain energy in “beam strips”. The method of Appendix F has proven 
very effective when applied to the 3 -node element [15,16], When applied to 
the 4-node element, it is too soft under some conditions, so an arbitrary 
reduction in the softening effect may be appropriate. Additional study of the 
application of the SABM to the 4-node element might lead to significant 
improvement. 


CONCLUDING REMARKS 

Some techniques for improvement of Mindlin shell elements for analysis of 
laminated composite aerospace structures, developed with support from 
several NASA contracts, are brought together in the Appendices of this 
report. These techniques have been applied to elements with 3, 4, 6, and 10 
nodes, and can be applied to other nodal configurations. Performance data for 
the MB4 4-node quadrilateral and the MB6 6-node triangular elements, 
which use these techniques, are also presented. 

The MB6 6-node triangle has uniformly excellent accuracy in all deformation 
modes, needs no options for different conditions, and performs much better 
than elements in some commercial programs in critical test cases. In regular 
rectangular meshes it is about equivalent to the 8-node quadrilateral with 2x2 
integration. In distorted meshes it is distinctly better. The triangular shape is 
an advantage in modeling. It should be universally adopted for stress analysis 
of laminated composite aerospace structures. 

The LU71 3-node [10, 16] and MB4 4-node elements have proven to be quite 
effective and robust, particularly in vibration analysis. 
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Figure 1. NAFEMS LE2 
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Figure 2. NAFEMS LE3 


SCORDELIS-LO ROOF [1,2J RESULTS 


Problem Statement 




Length: L = 50.0 Radius: R = 25.0 

Thickness: t =0.25 

Loading: uniform vertical gravity load of 90.0 per 
unit area 

Material Properties 

Isotropic: E = 4.32 x 10 8 , v = 0.0 

\Jr 

V s -! / 

\J/ 

r v 

V 

Boundary Conditions 

Supported on each end by rigid diaphragms, 

i.e. u = w = 0.0 on curved edges. 

r 


Target Output 

Vertical Displacement w at midside of free edge 
= 0.3024 


Figure 3(a) Scordelis-Lo Roof 
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Figure 3(b) Scordelis-Lo Roof Results 


Target Deflection = 0.3024 ft at mid-side of free edge 



Vertical Deflection at Mid-Side of Free Edge (ft) 

Nodes per Side 

MB6 

AIZ6 

AIZ8 

MB4SI 

3 

0.29178 

0.09344 

0.39091 

0.47804 

7 

0.30244 

0.22418 

0.30960 

0.31441 

13 

0.30495 

0.28884 

0.30364 

0.30572 

25 

0.30377 

0,29993 

0.30232 

0.30313 


Table 3. Scordelis-Lo Roof Results 
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Figure 4(b) Morley Skew Plate Results 


PLUNKETT’S VIBRATING WEDGE [7] RESULTS 



Figure 5. Plunkett’s Vibrating Wedge 


Objective 

Determine the first 12 out-of-plane modes of vibration for a cantilevered 
wedge section plate using a 3 x 6 mesh and consistent mass matrix. 

Material Properties 

Isotropic: E = 10E6, v = 0.3, p = 0.0002591 

Boundary Conditions 

All degrees of freedom are fixed at support. All u,v, and 9 Z displacements are 
fixed to eliminate in-plane vibration modes. 

Geometry 

a = 30, thickness at tip = 0.001, t = 0.96899451 


Results 

Target: Experimental Data (normalized) for mode shapes similar to 
Plunkett’s sketches. 


Element 

MB6 

AIZ6 

AIZ8 

MB4SI 

Mode 

Target 

Frequency 

Percent Error 

1 

2.47 

-2.8 

-2.4 

-3.0 

-3.1 

2 

10.6 

4.4 

6.9 

2.1 

-1.8 

3 

14.2 

1.4 

3.6 

1.1 

4.7 

4 

28.7 

-0.9 

7.7 

-1.7 

-5.6 

5 

34.4 

-0.6 

11.3 

-1.0 

-1.7 

6 

47.4 

-4.7 

13.9 

-4.5 

-9.5 

7 

52.5 

-3.4 

8.3 

-9.1 


8 

54.0 

-1.9 

18.1 

-3.6 


9 

63.5 

-9.1 

15.8 



10 

68.0 

-10.1 

27.6 




Table5. Plunkett’s Vibrating Wedge Results 
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Figure 6, NAFEMS T1 


Version 

Z Deflection 
At E (mm) 

% Error 

Bending Stress 
At E (MPa) 

% Error 

Interlaminar Shear 
Stress at D (MPa) 

% Error 

MB6 

-1.0545 

-0.5189 

683.14 

- 0. 1 1 1 1 

-4.00 

-2.4390 

MB 10 

-1.0544 

-0.5283 

684.60 

0.1024 

-4.01 

-2.1951 

AIZ6 

-1.0541 

-0.5566 

684.74 

0.1228 

-4.07 

-0.7317 

AIZ8 

-1.0544 

-0.5283 

684.23 

0.0483 

-4.10 

0.0000 

AIZ10 

-1.0542 

-0.5472 

684.28 

0.0556 

-4.10 

0.0000 


Table 7 Results for NAFEMS Laminated Strip 


NAFEMS 

LAMINATED STRIP 

Test No. 1 

draft 

Origin 

NAFEMS report 



Analysis type 

Ortbotropic 




Geometry 


y n 


0” fibre direction 
► 


L »* -i. 


15 

A. >« J 

r "T T *T* H 


lON/mm 



X A 

A 


A 

B 


~T 


AU dhoenshas m mm 


Loading 


Load line of 10N /mm at C (x-25, z-l) 

Boundary conditions 

One quarter model, simply supported at A (z=0) and 
reflective symmetry about x=25 and about y=5 

Material properties 

E, = 1.0E5 MPa, v l2 = 
G i2 = 3.0E3 MPa, v u 

0.4, Ej - 5.0E3 MPa, v, 2 = v 2l 

0J, G}j - 2.0E3 MPa Ej~ 

Element types 

Laminated beam, laminated plate, laminated Melt 
or stacked brick 

Meshes 


_ rn ' 

1 

1 — l - *^ 

1 j C, D, E 





Output 

Bending stress at E 
Interlaminar shear stress at D 
z deflection at E 

Target 483.9 MPa 

-4.1 MPa 
-1.04 mm 


Figure 7 NAFEMS Laminated strip 
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NAFEMS BENCHMARK TEST No. LE2 
For shells 

SYSTEM 

CASE 1 

CASE 2 


Nodes 
4 8 

Nodes 
4 8 

NASTRAN 

-153 

-10.5 

-17.0 

-6.0 


ANSYS 

-18.0 

+32.0 

-55.0 

+5.0 


GIFTS 



+2.9 


Membrane locking 

MELINA 

-98.2 

-9.8 

-2.3 

-5.0 


MELISSA 


-7.0 


-63 


MARC 

-4.2 

-12.0 

+10.8 

-6.9 


PAFEC 


+5.7 


-1.6 


ASAS 


-55.0 


+4.5 


LUSAS 

+4.0 

+3.1 

+7.6 

+7.6 


FINEL 


-10.0 


+2.0 


SUPERTAB 

-18.8 

-25.8 

■40.7 

■15.8 


BEPSAFE 


-2.2 


-8.2 

Reduced integration 

COSMOS/M 

-6.3 

4.8* 

-24.0 

+33* 


ABACUS 

0.0 

-0.3 

+1.0 

+0.8 


NISA 

+2.6 

-12.0 

+11.0 

-5.0 



NAFEMS BENCHMARK TEST No. LE3 
Hemisphere with pinch loading 

SYSTEM 

COARSE 

FINE 


Nodes 
4 8 

Nodes 
4 8 

NASTRAN 

+3.2 


+1.1 

+3.2 


ANSYS 


-19.4 

0.0 

+1.1 

Warp flag 

GIFTS 

49.6 


-7.2 



MELINA 


-5.4 


-8.1 

9-Node shell 

MELISSA 


-23.8 


-0.5 


MARC 

455.9 

-29.5 

-11.8 

-1.3 


PAFEC 


-113 


0.0 


ASAS 






LUSAS 

-93.5 

-9.7 

-81.6 

0.0 

Ftl integration 

FINEL 


-3.1 


-2.2 


SUPERTAB 

+1.6 

-31.9 

-0.5 

4.3 


BERSAFE 


-14.7 


-0.1 

Reduced integration 

COSMOS 

-11.3 

■15.0* 

-1.1 

-1.6* 


ABACUS 

+3.2 

4.3 

-1.1 

0.0 


NISA 

-9.5 

-11.4 

-173 

+7.6 



Table 8a. NAFEMS LE2 Curved Shell 
Patch Test 1989 Results 


Table 8b. NAFEMS LE3 Pinched Hemisphere 
Shell 1989 Results 
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ERROR SUMMARY TABLE - 6 AND 8 NODE ELEMENTS 

Percent Errors 

NASTRAN 
“B” Option 

-7.85 

(-12.68) 

© 

® « 
«s ^ 

-51.14 

+5.76 

+3.53 

990+ 


+3.72 

NASTRAN 
“B” Default 

30 

00 

+ 

+15.32 

(-53.3 to +15.32) 

09*e+ 

68‘0+ 

-2.83 

+9.55 

DIVERGES 


-17.75 

NASTRAN 
“A” Option 



+1.95 

+6.04 





NASTRAN 
“A” Default 

-3.41 

+85.33 
(-99.3 to +100) 

+16.43 

+15.40 

H-0.90 

® 

o 

o 

1 

+1.10 


AIZ8 

-10.97 

-6.20 

-25.16 

-0.23 

+2.38 

-0.03 

-35.20 

. . 

+2.28 

A1Z6 

i 

-12.08 

-0.46 

-91.98 

-59.24 

-25.87 

-0.82 

-18.40 

+1.89 

MB6 

-0.12 

-0.12 

— 1 

+0.49 

l 

-0.27 

+0.01 

+0.45 

+0.80 

91*1+ 


NAFEMS LE2 
Bending 

NAFEMS LE2 
Membrane 

NAFEMS LE3 
Coarse 

NAFEMS LE3 
Fine 

Scordelis-Lo Roof 
7 Nodes per side 

Scordelis-Lo Roof 
25 Nodes per side 

Morley Skew Plate 
211 d. o.f. 

NAFEMS LE5 


18 


Table 9. Error Summary for 6 and 8 Node Elements on Six Critical Test Cases 


Appendix A 

Curved Element Development 



The element x-y axis are in the plane defined by nodes 1-2-3, with origin at the centroid, 
x axis parallel to the 1-2 edge, and z axis perpendicular to the plane. The element matrices 
are developed in these axes. 

Unit vectors, e, and e 2 tangent to the mid-surface at any point are 


^ = >^ = X N ^y> 


e 2 


x ’7 H N i.n X ‘ 


[Al] 


Where x„ y h z/ are coordinates of nodes and N t are conventional interpolation functions. 
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k 

* h 



e x x e 2 
\e x xe 2 \ 



V l = V 2 xV 3 = 


m, 

L n x. 


[A2] 


The local axis x -y -z 'at any point are the directions of V { , V 2 , and V 3 , respectively. 


Rotations in element axis are 


X' 


' h 

k 

k' 

- 

a 


= 



”h 

fi 

X. 


_ n i 

n 2 

r h. 

0 


In a flat element 


[A3] 


V x = f, V 2 =j, a = 0 x and 0=0 y . 

With these definitions, the material property transformations can be done as in a flat 
element, provided that, the angle ^ which is input is the angle between the local x' axis 
and the fiber direction. 
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Curved Element Development 


* = £*<*,+! AT, 

<=l i=l ^ 

(=1 i=l ^ 

1=1 (=1 * 

Where n is the number of nodes in an element. 
Displacements in the element x-y-z axis are 

» = t N i u > + * ‘uPi i 

i=i & i=i 

v = Z + T^Z N i\r”hfl t +m u&] 

i=i ^ <=i 


[A4] 


[A5] 


w = Z + 7^Z + »uA ] 

i=i ^ i=i 


Where 

it n n 

Z^H 2 |« + / i ^ Z^if- /n 2i a + OT l</ ? ] Z^t _n 2/ a + ' I l<^ 

/=1 , (=1 , ;=1 


interpolates rotations of normal lines. This differs from Ahmad-Irons-Zienkiewicz in that 
interpolated nodal rotations are multiplied by the local thickness, /, scaled by 4" to produce 
a contribution to u, v, and w. This is significant only in variable thickness elements. 

Note that both rotations and displacements are defined at the middle surface, and the 
Jacobian J should be evaluated there. 

The transformation of rotations from the local axis a„ and /?, to the element axis Qy, 0 Z 
is 
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a i 


X 

"hi 

"u 

V 



Pi 

= 

k, 

"hi 

"2i 


= T0 

i t 

[A6] 

0_ 


lv 

"hi 

"y. 





The Jacobian J is 


/ = 


y*t z >t 

*», y > 7 z > 7 

■*»£■ z »f_ 


[A7] 


Where 


*»« = Z N i >s x i + Z N ‘ >s\ 9 k 




r=»l 


+i>,. Mj* 

/=! i=l ^ 


[A8] 


x K = lL N >\ t ht 


<- 1 


Etc. where n = number of nodes. 

In the new element versions stiffness matrix, the Jacobian / is always evaluated at the 
mid-surface, <£"= 0, so the terms containing % after differentiation vanish. 

Let (for 4-node case) (other cases similar) 


R u = 


R = 


k\ 

Ai 

0 

0 

0 

0 

0 

0 

0 


0 

0 

o' 

"7, 

0 

0 

o' 


0 

0 

0 

~k 2 

In 

0 

0 

0 

0 


0 

0 

0 

0 

T 2 

0 

0 


[A9] 

0 

0 

0 

0 

0 

0 

~^23 

In 

0 


0 

0 

0 

0 

0 

r 3 

0 


_ 0 

0 

0 

0 

0 

0 

0 

0 

0 



'u 

0 

0 

0 

0 

T*. 


-«u 

"u 

0 

0 

0 

0 

0 

0 


0 

0 


0 

o' 

X 

0 

0 

o' 


0 

0 

0 

- n 22 

n n 

0 

0 

0 


0 

0 


0 

0 

0 

T 2 

0 

0 


0 

0 

0 

0 

0 

0 

-"23 

"n 


0 

0 


0 

0 

0 

0 

T 3 

0 


0 

0 

0 

0 

0 

0 

0 

0 


0 

~"2A 

”14 

0 

0 

0 

0 

I 
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'-"*21 

"hi 

0 

0 

0 

0 

0 

0 

0 

0 

0 

O' 

X 

0 

0 

0" 

0 

0 

0 

-"hi 

"hi 

0 

0 

0 

0 

0 

0 

0 

0 

T 2 

0 

0 

0 

0 

0 

0 

0 

0 

~"hi 

" \3 

0 

0 

0 

0 

0 

0 

T3 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

~"h4 

"h4 

0 

0 

0 

0 

T 4„ 


Therefore, the displacements within our element can be expressed as a function of the 
nodal displacements by 




N 

0 

0 

1 

es a 

1 < 

u 

u 

V 

_ 

0 

N 

0 

l 

V 

w 


0 

0 

1 V 

2 

w 

e 


[A10] 


Where, for the 4-node case, 
u r = [u x u, u, u 4 ] 

^2 ^3 n] 

[All] 

w 2 W 3 W 4 ] 


® - \Oxl Qy\ e zi @x 2 @y 2 @z 2 ^xi & y 3 @z 3 @ x 4 & y 4 ^rt] 

To get the strains we differentiate u, v, and w, for example 



r 1 1 

u 

V 


r 1 i 

u 

V 


N,< 0 0 

w 

0 

+ 

0 0 0 -Ct, f NR u 

w 

0 


[A12] 


Terms involving the derivatives of t such as the t, § have been neglected since they give 
strains due to rigid body motion. These derivatives vanish when the element thickness is 
constant. This formulation is used to avoid singularity due to very thin edges. The 
Ahmad-Irons-Zienkiewicz formulation may be a bit better for thick elements. 
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The derivatives of the displacement functions may be computed at any point as 



N,, 

0 

0 






0 

0 

jCi".,*. 

1 




0 

0 

0 

2 tNR * 


u >< 


0 


0 


u 

v 

% 

= 

0 


0 

V 

W 

V V 


0 

0 

0 

\i‘N, e R. 

9 



0 

0 








1 



0 

0 


2 & N >r, K 



0 

0 

0 

\' NR - . 



[A 13] 


Note that this matrix can be divided into a part including 4" and a part not including 4" 




J ■* 

0 

0 


0 

/ 1 

0 


0 

0 




w 


’C 


[A 14] 


Matrix H above transforms derivatives to strains in the element {x-y-z) axis (Cook, 
360) and matrix F e transforms strains to local (x -y -z ) axis (Cook, p. 212). 

Combining these equations, we can write 


U 

V 

w 


p- 


[A15] 


9 
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Where strains are in the local x-y'-z 'axis. 

The stiffness matrix (integrated in local axis) is then 


K = \B T D'B[DetJ]d^d7jd<; 


[A 16] 


Only the stretching and bending parts of the stiffness matrix are computed from this 
expression in our elements. The out-of-plane shearing stiffness is computed from a 
different process. 

In these equations, the Jacobian, J, is evaluated at the middle surface. 

Note that explicit integration through the thickness is possible. 

JT = J [Bl + (IB[]D'[B, l„ + OBjDetJ}i&nd( 

= j B’ 0 D’B„[Detjy&iidt;+\ B![f ! t 2 D']B,[DelS}i#/r/d( [A17] 

+ 1 [Si [<?!>’]«, + B,‘ [0IT]B'lDet 


B 0 contains /, but in the form above, l is factored out of B\. 

J><= 




= D[ 


,c 


2l^> 


- 1-1 


r 2 

+ D'2- 

2 -* 


NLAl 1 

-Z 

7=1 Z 


[A18] 


In a sandwich structure with constant-thickness skins and variable-thickness core, the 
integrals 

j'_D'd(, andjVfVf 
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would have different values at each integration point. 
To complete the integration of 


x-j #[/', D'd^BJJet[j]d^ 

+ J D' [A1S 

+ J [« 0 r /[jV «/<■]«, + D'&c]B,]Det[j)d& n 

does not require multiple points in the through-thickness direction. The computing time 
could be considerably reduced, compared to using two points per layer in the through- 
thickness direction. 
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Appendix B 

Modifications to 

In curved quadratic elements, bending causes spurious mid-surface strains, except at the 2 
x 2 Gauss points. These spurious strains cause serious errors in fully integrated Ahmad- 
Irons-Zienkiewicz elements when the “rise” (deviation from flatness) of the element is 
more than about 1/5 of its thickness. 

It has been observed that, in constant-curvature bending of quadratic elements, the 
average mid-surface strains is (correctly) zero. Thus, the spurious mid-surface strain 
which causes excessive bending stiffness can be avoided by replacing B 0 with 


where m is the number of integration points over £and rj and ff* are the integration 
weights. 

Unfortunately, this eliminates the gradient in mid-surface strains which is needed for 
accurate solution to some problems. It also introduces mechanisms. 


flat 6-node triangle with straight sides. No mechanisms are introduced by this process. 
Strain-displacement matrices are computed for the 3 comer subregions as flat 3- 
dimensional constant strain triangles. 


m 




[Bl] 



/ \ 


In the 6-node triangular element, the 
gradient can be restored by using 
triangular subregions. Using 3-point 
integration for the 6-node triangle, the 
integration points coincide with centroids 
of the 3 comer subregions.' Denoting the 
strain-displacement matrices of the three 
subregions by B', B* 2 , Bl , the membrane 
strain gradient is restored by using: 



It can be shown that this gives exactly 
correct membrane strains for all constant 
strain states and all linear strain states in a 
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Appendix C 

Shear Perpendicular to Surface in Flat Elements 


In Mindlin plate theory, shear perpendicular to the surface 


Y 


yz 



is given by 


dw 


(C.l) 


Sign conventions for displacements and shear strains are shown in Figure C.l. 

It is evident from equations C. 1 that the function representing w should be one order 
higher than the functions representing©* and 0 y . Some conventional element nodal 

configurations and the terms appearing in interpolation functions are shown in Figure C.2 The 
solid lines bound terms in- the conventional interpolation functions used for©* and 0 V , and the 


dashed lines bound terms required in the “unconventional” interpolation functions used for w by 
Family 1 elements. 


All elements considered here have three degrees of freedom at each node: displacement 
w, rotation 0* , and rotation 0 V . Rotations are always represented by conventional 

interpolation functions which are continuous on and across element boundaries. Bending strain 
energy is computed from rotation interpolation functions in the conventional manner. 


Element Family 1 

Following Utku (1) we choose 
w = w' + w' 


(C.2) 


in which w' represents bending or “Kirchhoffian” deflection and w* represents shearing 
deflection. 


dw _ dw' dw" 

3 3 3 

dw _ dw' dw 

dy cty dy 


(C.3) 


When shear deformation is zero we see from Figure C. 1 that 
dw _ ^ _ dw' 

~3~ " y ~~3 

dw _ q _ dw' 
dy " x dy 

Using C.3 


(C.3a) 
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(c) Shear deformation without bending 
Figure C. 1 Sign Conventions and Shear Deformation Relations 
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Nodal 

Configuration 



Terms Included 
in Interpolation 


1 


* y 

x 2 xy y 2 
x’ x’y y 2 x y ’ 



1 


x y 

x 2 xy y 1 


x‘ x’y y’ x y 


) 


x 4 x’y x’y 2 y’x y 4 


Number 

Equations 

Available 


Number 
Unknowns 
Family | 


6 


5 


12 


9 



1 

x y 

x 2 xy y 2 
x 2 x’y y’x y s 

x 4 xy x’y 2 xy 1 y 4 
x’ x 4 y x’y 2 x’y 2 xy 4 y 1 
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14 




1 

x y 

\x’ xy y 2 / 
x \j’y y’x/ y 
x\x‘y x’y’ xy/y 4 


X 5 xV xY xY xy 4 y s 


16 


12 



Figure C.2 Nodal Configurations and Interpolation Terms 
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(C.3.b) 


dw* dw 
3c 3c 


+ ®y=Y* 


dw’ dw „ 

’’ r ’ 


In Figure C. 1 it is shown that 


3w 


Y a w hen — = 0 


3c 


dw 


r =-~ when 0=0 
3c 

so the choice of functions in Eq. C.2 is acceptable. 

The method of Element Family 1 will be illustrated by application to the four-node 
element. 

The conventional, continuous functions used for 0 are 
0, = <% + a x x + ay + a^xy 

Q y = b 0 + b x x w-by* b r xy 

The function for w ' is chosen to contain terms one order higher. 

w’ = e 0 + c,x + c 2 y + c 3 xy + c 4 : d + c 5 y 2 + c 6 x 2 y + c 7 y 2 x 

An immediate problem arises, since differentiation of w' does not yield exactly the 
desired expressions for 0^ and © v 

dw’ 

0* = = c 2 + CjX + 2 cy + c 6 x 2 + 2c 7 yx 

dy 

The coefficients can be identified with those in Eq. C.4 except for term c 6 x 2 . 


(C.4) 

(C.5) 


(C.6) 


-e 

' dx 


= c, + cy + 2 c 4 x + 2 c 6 xy + c 7 y l 


Again, coefficients can be identified with those in Eq. C.4, except for term c 7 y 2 . 

Another conflict which arises is that c 3 corresponds to both a, and b 2 in Eq. C.4. Utku 

resolved this problem in his three-node element by using an average value. 

Numerous attempts have failed to find a way to eliminate these conflicts. However, it 
appears that they do not prevent reaching the goal of representing all constant strain states 
exactly. 

For pure bending in the y-z plane 
c^w’ 

— r- = 2 c 5 + 2 c 7 x = constant 


(C.7) 


dy 2 

which requires that c 7 -0 for this pure bending state. For pure bending in the x-z plane 
d*w' 

— r- = 2c, + Icy = constant 
3c 

which requires that c 6 = 0 for this pure bending state. 

For pure twist 


(C.8a) 


(C.8b) 
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c^w' 


= c 3 + 2c 6 jc + 2c n y 


dxdy 

which requires that c 6 and c 7 = 0. 

To compute the values of the c j , we require that Eq. C.6 and Eq. C.7 be satisfied at the 
nodal points. The bar above a symbol designates nodal values of the variable. 


(C.9) 
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Denoting the 8x7 matrix above as A 

© = Ac 

A r © = A r Ac 

c = [A r A] _1 A r 0 

This is a least-squares solution for c. 

W' = C„ + |X 


(C.10) 


y xy x 2 y 2 x 2 y xy 2 ][A r A] A 7 © 

Note that additional terms could be added to w' , and Eq. C. 10 can be solved if Matrix A has at 
least as many rows as it has columns. 

Interpolating w* by the same function as 0, and 0 V . 

w* = d 0 + d x x + dyy + d 2 xy 

w = w* + w' 

Combining c 0 with d 0 

w = [l x y x^]d+[jc y xy x 2 y 2 x 2 y jcy 2 | [A r A] 1 A r © 


(C.ll) 

(C.12) 


(C. 13) 
(C.14) 
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w = Bd + C[A r A]" l A r © 

B r W = B r Bd + B r C [A r A] _I A r 0 
d = [b'bJVw - [B t B]‘b 7 C [a’aJ'V® 

w*=[l jc y xy]d-c 0 

w = [l x y jy]d+[.r y xy x 2 y 2 x 2 y xy 2 jc 
w = [ 1 * y ^y]|[B r B]' l B 7 W-[B r B]' ! B r C[A r A]' I A r ©| 


(C.15) 


(C.16) 


+[x y xy x 2 y 2 x 2 y xy 2 ] [A r A]~ A r © 
Shear strain can be computed using Eq C.3b 

a*' 


r s = 
r« = 


dx 


• = [0 1 0 y] d 

= [0 0 1 x]d 


& ' ' (C. 17) 

An alternate way of expressing w* is, using the conventional isoparametric interpolation 


functions 


w* = N w’ 

where W* represents values of w* at the nodal points. Equation C. 14 can be re-written as 
w = w* +[jc y xy x 2 y 2 x 2 y xy 2 ][A r A] 'a 7 © 

W = w* + c [a 7 a] V® 


(C. 13a) 
(C.14a) 


Then 


and 

w* = w- C [A r A] ‘A 7 © 
w ’ = nJw-C [A r a] ’ A 7 ©] 


(C. 15a) 


(C.l6a) 
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d*v ek N 


|w-C[A r A] A 7 ©! 
|w-C[A r Aj ’a 7 ©] 


Sc Sc 
dw ' fit 


* dy dy 


(C.17a) 


Equation C.17a provides a simpler way of calculating shear strains. The expression 
inside the brackets contains the nodal deformations but otherwise consists of constants. 


Element Family 2 

Calculation of the coefficients in the expression for w' is done by writing equations 
which state that at nodal points 


Thus, each nodal point provides two equations. In Figure C.2, the last two columns show the 
number of equations available for calculation of coefficients and the number required for Family 
1 elements. 

For the four-node element, eight equations are available, and there are seven coefficients 
to be calculated for the Family 1 element. One additional term could be added. Since symmetry 
is required, the only choice is to add the x 2 y 2 term. 

Some of the other elements in Figure C.2 offer several choices of additional terms. The 
six-node triangle allows addition of either one or three terms. The eight-node quadrilateral 
allows addition of either two or four terms, but there are three possible feasible patterns. 

Calculation of element matrices for elements of Family 2 proceeds exactly as for 
elements of Family 1. There are more terms in w' defined in Eq. C.5; more columns in Matrix 
A, defined in Eq. C. 10; and more unknown coefficients in vector c, but the calculation process is 
otherwise identical. 


Out-of-Plane Shear Using Weighted Least Squares 

In the generalized Utku procedure for generating the out-of-plane shear stiffness submatrix, the 
nodal rotations, 0 , are related to parameters c from the function describing the “Kirchhoffian 
deflection”, w ' by 


Usually, matrix A has more rows than columns, and a least squares solution for c is required. 
This solution can be weighted by introducing a diagonal matrix of weights, W, For example, 
side nodes of a 6-node triangle could be weighted differently from comer nodes. 



(Cl 8) 


0 - Ac 


[C. 19] 
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[C.20] 


W® = WAc 
A r W® = A r WAc 

c = [a t wa]' 1 a t w® 

This weighting can be easily introduced by replacing A T by A T W. 

For example, to weight comer nodes differently from side nodes in the 6-node triangle, with 
rotations order as 


Simply multiply the first 6 columns of A T by the comer node weight coefficient, which leaves 
W t = 1 for side nodes. An extensive numerical study indicates that 0.07 is a good value for comer 
node weights in 6-node elements. 

This treatment of shear strains is strictly valid only for flat elements, but works very well in 
curved elements if the included angle between outward normals is not too large, as demonstrated 
by test cases. In the Scordelis-Lo roof it works very well with only two elements, each spanning 
40 degrees. However, most performance data is from standard tests where the elements only 
span about 20 degrees. 



[C.22] 
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Appendix D 

Strain Calculation at Nodes 


Strains are initially calculated at integration points, but values are desired at other 
locations, especially at nodal points. Strains are derivatives of displacements, and should 
be interpolated by functions one order lower than displacements. 

In the case of the 6-node triangle, displacements vary quadratically, and linear 
functions are used to interpolate strains. Using 3 integration points, the relation between 
integration point strains and strains at the 3 comer nodes is 
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Strains at any point can be obtained as 
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or, transposing 
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This relation can be used to compute strains at any point, including the 6 nodal points. 
The same process is used for all components of strain; membrane, bending and shear. 


Interlaminar Shear Stresses 


Interlaminar shear stresses can be obtained from integration of the equilibrium equations 
da, dr. 


dx 




dy dz 


dr„ da„ dr 
v - + 


dx dy 


'+— *- = 0 


\D6] 


dz 


dr„ dr da 
dx dy dz 

We compute in the local axis 


= 


r 

r ( ar xy 3 

J -//2 

t ax’ ay' J 

r 

f 5<r . dr . , 3 

J -//2 

^ ay' dx' ) 


dz' 


[D7] 


G x\x‘ 


S x\x' 




1 

w. 

V 

1 

G y\*' 

= /> 

*** 

and 

<Xyy 

= Z) 

Sy.y 

X r . , 

JC > _ 


Y x'y\x\ 


r , , , 

L *y,y j 


Y x'y ,y_ 


PM] 


In Mindlin plates and shells, the strains must vary linearly through the total thickness, and 
stresses must vary linearly within a layer where material properties are constant. 

* 



top of layer k 
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f yy ] 1 
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Where a x , x , and r xy , are average or mid-thickness values in layer / and t, is thickness 
of layer i. 

Similarly, 



top of layer k 


= + r xy *i 

i-1 i 


PIO] 
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The ANSYS theoretical manual notes that this summation which starts with r = 0 
and r x , : , = 0 at the bottom may give t x ,,, * 0 and or r y , z , * 0 at the top due to inexactness 
in the calculation. ANSYS has a procedure for distribution this error to make r x , z , = 0 at 
the top, which is also used in our programs. 

In the 6-node triangle, which has 3 integration points, strains needed are generated as 
shown below, and then multiplied by material properties as in Eq [D8] to produce the 
stresses needed in Eq [D9] and [DIO], 
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In the 6-node triangle, linear interpolation functions for strains, N t , are equal to area 
coordinates <§, and derivatives in Eq. [D1 1] are 


^1,*' =C*' = AC +/ ”lC 

N 2,S=$2,x'=h r l,x+ m x r l.y 

N 3 ,x' = = "AC " ™iC - A >7,, " ”W. y 

AA y = C/ = AC+^C 
Cy = C = A^ +m 2^ 

N x/ = Cy = ~4C - ”hly - kv, x - "hJl, y 


[D12] 


Appendix E 

Out-of-Plane Rotation Stiffness 

Mindlin plate and shell elements usually have no rotational stiffness associated with the 
component of the nodal rotation vector that is perpendicular to the element surface. (An 
exception to this is elements with Allman rotations) 

All of our plate/shell element matrices are formulated in element axis with 6 degrees of 
freedom per node. Two procedures are incorporated which can provide the needed out- 
of-plane rotational stiffness. Each has a coefficient that can be an input item in the 
program, to modify the “basic” stiffnesses. The combination of these procedures avoids 
singularity in flat regions without fixing any degrees of freedom. 

Zienkiewicz procedure adds a stiffness to the diagonal of each 0 Z degree of freedom, and 
subtracts l/(n-l) (n = number of nodes) of that stiffness from each off-diagonal 0 Z degree 
of freedom in the same row and column. This retains symmetry, equilibrium, and 
freedom from forces under rigid body motion. Our “basic” stiffness is the average of the 
nodal in-plane rotational stiffnesses from the bending sub-matrix. 

The Kanok Nukulchai procedure is a penalty method which creates an artificial quasi- 
strain energy, U,. 



[HI] 


The quasi-strain, e h can be written as 


e, 



[E2] 


Where 0 . local rotation perpendicular to the surface interpolated from 

nodal values 
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elasticity rotation 

In-surface shear modulus 
thickness 

Coefficient for magnitude of stiffness 
vector of nodal displacements 


Then, the artificial stiffness matrix is 

K, = * B?B,dA *c,G, y f'W l [B[B,]DelJ 

1 = 1 

Where n = number of integration points and w, = weights. 
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The determinant of the 3D Jacobian used in our programs is evaluated at the mid-surface. 
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For a 4-node element 
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An extensive numerical study indicates that 1 percent of the average of the in-plane 
rotational stiffnesses from the bending submatrix is a good value for the Zienkiewicz 
procedure, and that 0.05 is a good value for the coefficient in the Kanok-Nukulchai 
procedure. 


Appendix F 

Simulated Antisymmetric Bending Mode 

Mindlin elements with 3 or 4 nodes are not capable of anti-symmetric bending, i.e. 
correct response to linearly varying bending moment. The simulated anti-symmetric 
bending mode (22, 23, 24) substitutes shear deflection for bending deflection and 
preserves the correct total strain energy in “beam strips”. 

For one layer, a reduced shear modulus, G^ , for a strip in the local x direction is 


i i ,Kf 

g; g m E x h 2 


[Fl] 


Where L x is the length of the strip, E x is the bending modulus, and h is the thickness 
(depth) of the strip, and a is an effective length factor. 

For symmetric orthotropic layers 
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For general orthotropic layers 
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The notation for Q is used as defined by R. M. Jones. 
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In the local x axis, C55 = G s z and Q u = E x . 


In the four node 
element, 4 strips are 
used. The effective 
length factors, a, are 
0.707 for the diagonal 
strips and 1 .0 for 
strips connecting 
midsides. 
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Values of j C* 5 on the left hand side are calculated for each “beam strip”, which requires 

appropriate material property transformations. The equation above is then solved for by 
least squares to yield 

jc^dz, J Cl $ dz, and J C' i} dz 

in the element axis. 

If the shearing stiffness sub-matrix is integrated through the thickness before integration 
over the area, these integrals are used directly. Otherwise, ratios of reduced integrals to 
true integrals are computed and these ratios applied to out-of-plane shear moduli. 


In the triangular element, strips are taken parallel to each side and the theoretical 
effective lengths should be a- 0.707. 
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This procedure has worked very well in 3-node triangular elements. In 4-node elements it 
is very good in some deformation modes, such as in a cantilever beam. However, in a 
centrally loaded square plate modeled by 4-node elements it is too soft, so the 
coefficients need to be reduced from the theoretical values derived above. More study of 
this is needed. 
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APPENDIX G 


6-Node Flat Shear Panel (Using 10 Node Shape Functions) 

Out-of-plane shear strains can be calculated entirely in curvilinear coordinates, as shown here. Numerical 
results from this process are identical to those obtained by using a complete cubic polynomial in the 
procedure of Appendix C. While the development in this Appendix is strictly valid only for a flat shear 
panel, it is an important step towards development of a curved shear panel. The interpolation functions for 
a 10 -node triangle are used as a convenient representation of a cubic polynomial, and the parameters, C/, 
just happen to correspond to deflections at the nodes of a 10-node triangle 

The total displacement w from bending and shear is expressed as 

w = w* + w + [Gl] 

Where w>* represents shear displacement and w + is bending. 

Note that 

W + = N,C, [G2] 



Figure G-l 10-Node Triangle 

The functions, N j, are 
For the comers (/= 1,2,3) 

N,=j£(3f,-lX3£-2) 

For the mid-side nodes and node 10 
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N 10 - 276ft* 

N 4 =|ftft (3ft -1) 

N, = |ftft( 3ft - 1) 

N s =|ftft(3ft-1) 

N,=|ftft(3ft-1) 

N,=|ftft(3ft-1) 

N, =|ftft(3ft-l) 

We require H' + = N,C ( as a function of ft , ft , ft and c,. TWs allows us to compute 
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Finally, the derivatives of the functions N with respect to are 
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N,.. = i[27£- 18f,+2] 

N« = 0 

N l{ = -i[27^-18ft+2] 

N„ = |[««,ft-*] [ 08 ] 

N„-f[3g-&] 

N«— Ipg-fe] 

N,.. = |[f 3 (3f 3 - 1) - f,(3f 3 - 1) - 3f,f 3 ] 
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d % drj d% dr] 
We need — , , ~ 

ox ax oy ay 

Use 6-Node Shape Functions 
x = NjXj 

j = 1 to 6 
From Cook 
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Which allows us to solve for c as 


Note: we generate matrix A in 2-row blocks 
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The matrix A r A in equation G16 is initially singular. The singularity can be removed by augmenting 
matrix A with an equation that assigns a zero value to w* at the centroid of the element. After some 
manipulation, the only change is to add one (1) loA T A 10-10 which removes the singularity and allows 


inversion. 


The terms 
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in equation G17 can be obtained from the T in equation G13 defined earlier. 
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Which can be written more compactly as 

W = m m +X[A T AY X A T \e\ [G20] 


This can be rewritten as 

w* = w-X[A r A]- l A r [0] [G21] 


Finally, from basic mechanics we may write 
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Applying to the 6-Node flat shear panel yields 
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